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Stochastic dynamics govern many important processes in cellular biology, and an underlying 
theoretical approach describing these dynamics is desirable to address a wealth of questions in biology 
and medicine. Mathematical tools exist for treating several important examples of these stochastic 
processes, most notably gene expression, and random partitioning at single cell divisions or after a 
steady state has been reached. Comparatively little work exists exploring different and specific ways 
that repeated cell divisions can lead to stochastic inheritance of unequilibrated cellular populations. 
Here we introduce a mathematical formalism to describe cellular agents that are subject to random 
creation, replication, and/or degradation, and are inherited according to a range of random dynamics 
at cell divisions. We obtain closed-form generating functions describing systems at any time after 
any number of cell divisions for binomial partitioning and divisions provoking a deterministic or 
random, subtractive or additive change in copy number, and show that these solutions agree exactly 
with stochastic simulation. We apply this general formalism to several example problems involving 
the dynamics of mitochondrial DNA (mtDNA) during development and organismal lifetimes. 


I. INTRODUCTION 


Stochastic dynamics underlie a multitude of processes in cellular biology [THS]. Understanding the sources of this randomness 
C^ithin and between cells is a central current challenge in quantitative biology [7]. Noise has been found to affect processes 
mcluding stem cell fate decisions [HE], bet-hedging in bacterial phenotypes [ini E] , cancer development [12] , and responses to 
^i^poptosis-inducing factors HKH], illustrating the fact that a theoretical understanding of stochastic cellular biology is of great 
^importance in medical and biological problems. 

Partitioning of cellular components at cell divisions provides a considerable source of stochasticity in cell biology [TS| . Huh 
^nd Paulsson have shown that uneven segregration of cellular constituents at mitosis can contribute signihcantly to cell-to- 
. il^ell differences in levels of cellular components and proteins in a population, focusing on stochasticity in protein inheritance 
^^etween sister cells m- In addition to variability in protein levels, there is evidence that random partitioning of mitochondria 
cell divisions can lead to substantial extrinsic variability in the physical and chemical attributes, and behavioural phenomena 
i including differentiation propensity, in a population of cells [HI- 

Historically, mathematical modelling, including the use of birth-and-death processes, have provided a theoretical foundation 
T^ith which to describe phenomena in stochastic biology m- Early work on the problem of the stochastic evolution of cellular 
J^onstituents in a population of dividing cells was performed in the context of protein levels in bacterial cells by Berg [19j and 
.i;;^igney j^D], who made analytic progress with birth-and-death models coupled to cell division in the case of a steady-state 

,_population of cells with constant birth and death rates. A famous example of stochastic analysis of cellular systems is the 

y Qy iiria-neibriick treatment of population statistics of bacterial mutations [U], which has been addressed by several analyses 
Ohcluding generating function approaches [22]. Matrix equations, with operators corresponding to the processes of birth, death, 

,_partitioning, have also been used to obtain numerical results on stochastic effects in populations of dividing cells [231 [33] ; 

f*~a.nd a framework of nonequilibrium statistical mechanics has been used to derive general properties of protein content of dividing 
l/7)ells [23]. Stochastic models examining the behaviour of a continuous variable (for example, protein concentrations) have also 
’“^een widely used in cellular biology |26j . 

^ Several stochastic models have been formulated for systems involving compartmentalised elements which replicate and are 
• 'partitioned as compartments divide. An early example of this approach is due to Dowman [27], which considers mean copy 
l/^umber and extinction probability of cellular elements for specific birth, death, and partitioning dynamics. Other applications 
^re in the study of growing intestinal crypts, in which a collection of crypts, each containing a replicating quantity of stem cells, 
grow and divide according to a branching process |28| , and in the stochastic corrector model, whereby a population of replicators 
with given concentrations of constituents divide and propagate |29j . More recently, Swain et al. |30| have derived analytic results 
describing stochastic gene expression in dividing cells after a steady limit cycle has been reached. Quantitative results, in terms of 
integrals over kinetic rate dynamics, have been obtained for stochastic gene expression where cellular components are binomially 
partitioned at cell divisions m- The variance resulting from more general partitioning dynamics of cellular components has also 
been addressed for single cell divisions [I31[IS]. 

In this article, we focus on birth-immigration-death (BID) dynamics rather than the well-studied model dynamics of stochastic 
gene expression. With this framework we attempt to provide a model for the stochastic dynamics of cellular components other 
than gene products; we particularly focus on mitochondrial DNA (mtDNA) in several examples. Populations of hundreds or 
thousands of mtDNA molecules are typically present in eukaryotic cells, replicating and degrading somewhat independently of the 
cell cycle. MtDNA encodes vital bioenergetic machinery, making it an important target for stochastic analysis. We will consider a 
variety of partitioning regimes and, where possible, an arbitrary number of cell divisions, and aim to derive closed-form generating 
functions describing the dynamics of our model cellular components. In so doing, we avoid assumptions about equilibrium 
behaviour and restrictions to lower-order moments of copy number distributions, aiming to produce a non-equilibrium theory to 
describe stochastic dynamics in full distributional detail, at arbitrary times during any cell cycle. Our specific consideration of 
mtDNA dynamics provides a theoretical framework with which a class of models, often analysed numerically through simulation, 
can be descibed analytically. 


2 


A. Dynamics between cell divisions 


Random 

processes: 


0 


• Immigration m ^ m+1, rate a 

• • Birth iTi ^ m+1, rate Am 

0 Death m ^ m-1, rate vm 


B. Partitioning at cell divisions 

Exampie 1. Binomiai partitioning (fission-like division) 

(e.g. m = 6) 


• • 


• • • 

* • • 


• • 
• • 


• • 


• • • 


Exampie 2. Subtractive partitioning (budding) 


(e.g. m = 5) 


• • • 


• • 


• • 


(e.g. m = 11) 


' • 

• • 


• 

• • 


(e.g. m = 8) 


FIG. 1: Illustration of birth-immigration-death dynamics with random partitioning. We will derive expressions for the statistics of copy 
number m of cellular agents over successive cell generations, separated by divisions. (A) Between cell divisions, agents may be produced, 
replicated, or degraded; each is a Poissonian event. The copy number m of agents in a cell is a random variable that changes with these 
dynamics. (B) At cell divisions, the copy number of agents changes according to a different type of random event. Two possibilities are 
illustrated here: the binomial partitioning of agents into two daughter cells (one of which will be tracked in the next generation), and the 
loss of a random number of agents to a smaller bud (the larger remaining cell will be tracked in the next generation). 


In the first section we introduce our formalism and derive generating functions for BID dynamics with binomial partitioning 
and the inheritance of a deterministically or randomly reduced or increased complement of the parent cell’s population. We next 
illustrate the exact agreement between our theory and stochastic simulation, and investigate several example biological questions, 
regarding the dynamics of mtDNA during a copy number bottleneck and over organismal lifetimes involving many cell divisions. 
We further harness the power afforded by a generating function approach to explore extinction probabilities in these systems, 
and also perform our analysis for a given number of cell divisions in which cellular components are deterministically partitioned, 
and randomly partitioned in clusters. Finally, we discuss the implications of our mathematical formalism for approaches to 
stochastic biology. 


II. MODEL AND ANALYSIS 

In this section we detail the approach we take to obtain generating functions describing the stochastic dynamics of cellular 
agents subject to birth-immigration-death dynamics and stochastic partitioning at cell divisions. We first illustrate how the 
generating function describing agent dynamics within a cell cycle is derived. We next consider how this solution may be 
extended over cell divisions, using an assumption (later shown to be true) about the functional form of the expression describing 
the inheritance of agents at cell divisions. We then show that this extended solution gives rise to an overall generating function 
containing factors that are the solutions to recursion relations, where each recursive step corresponds to a cell division. We 
obtain solutions for these recursion relations, thus yielding general generating functions for dynamics over arbitrarily many cell 
divisions. Finally, we validate our early assumption for several important specific cases of inheritance dynamics and derive the 
specific solutions in these case. 


A. Agents within a cell cycle 

We consider birth-immigration-death (BID) dynamics, where agents are created (0 —•) with rate a, replicate (• —••) with 
rate A, and are removed (• —)■ 0) with rate v (see Fig. These processes are referred to as immigration, birth, and death 
respectively, and are assumed to be Poissonian, with non-time-varying rates (though see Results). We will later consider setting 
some of these parameters to zero, as special cases of the overall BID dynamics. We note that some literature refers to our 
‘immigration’ term (producing agents from nothing) as ‘birth’, using the term ‘replication’ to describe the production of agents 
from existing agents; the specihc symbols used to denote these rates also vary. For compatibility with a wide body of literature 
we adopt the nomenclature above. 

The dynamics of populations under BID dynamics are given by the corresponding master equation, describing the time 
evolution of the probability P(m,t) that the system contains m agents at time t: 


dP{m, t) 
dt 


aP{m — l,t) + v{m -\- l)P{m -I- 1, t) -I- A(m — l)P(m — 1, t) — (a -I- vm + Xm)P{m, t), 


with initial condition 


( 1 ) 


Pijyi-i 0 ) ^mmQ 


( 2 ) 
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We are concerned with the generating function G{z, t) = J2m=o t) for the distribution of cellular components in a gen¬ 

eral set of conditions. Once this generating function has been calculated, then, by its definition, all information about the distri¬ 
bution P{m,t) can be obtained through taking its derivatives, for example, E(m) = V(m) = ^ ~ (^)^) ’ 

P{m) = generating function corresponding to Eqns. 



dG{z, t) 
~~dt 

G{z,0) 


a(z - l)G(z,t) + (i^(l - z) + 

oz 

^ 1 


the solution to which is well known (see Appendix): 


/ v-X f ue'-^-'''>\z - 1) - Xz + v\ 

—1) — Xz + v) \AeO-‘^)*( 2 ; — — Xz + vj 

'■ -V--V-^ 

iGP gG,t) 

^ az,t){giz,t)ro 


(3) 

(4) 

(5) 

( 6 ) 


where the final line, using the underbraced substitutions in 


Eqn. casts the solution in a form that will be useful later. 


B. Partitioning of agents at cell divisions 

We now consider a series of cell divisions, linked by quiescent periods governed by the within-cell-cycle dynamics above. Time 
is accounted for throughout this model in the following manner. A full cell cycle is assumed to take time t, after which a division 
occurs. A variable t measures the elapsed time since the most recent cell division. As divisions occur every interval r, t is always 
less than or equal to r. If n divisions have occurred, the total elapsed time since the initial state is nr -I- t. We will here assume 
that r is constant, showing later that this assumption can be relaxed when our approach is extended to deal with many different 
dynamic phases. 

We will generally write Pi(rn, t|mo) for the probability of observing m agents at a time t after the ith cell division, given initial 
condition toq (at the start of the cell cycle before the first cell division). Hence, Po{m,t\mo) is the probability of observing copy 
number m at time t, given that zero cell divisions have occurred; and Pi-i(m,T\mo) is the probability of observing m agents at 
a time r after the {i — l)th cell division (i.e. immediately before the ith cell division). 

Consider the Ah cell division in a series of divisions. Throughout this article, we will use subscript a to denote ‘after’ and 
subscript b to denote ‘before’ a cell division:, thus, we write mi^b for copy number before the division and mi^a for copy number 
afterwards. We assume that mi^b > always. We write Ps{a\b) for the probability of observing a agents after a cell division, 
given b agents before that division. The generating function at time t after the cell division is then given by 


Gi{z,t) 

OO 'ITT'iyb 

= ^ ^ ^ z"‘Po{m,t\mi^a)Ps{'rni^a\mi^b)Pi-iimi^b,T\mo) (7) 

m mi^b—Orrii^a—O 
OO '>Pbi,b 

= E E G{z, t\m^^a)Ps{mi^a\mi,b)P^-lirni^b, T\mo) (8) 

rrii^b—O rrii^a—O 
OO 

= E E £,{z, t)g{z, Ps{mi^a\mi^b)Pi-i{mi^b, T\mo) (9) 

rrii^b—O rrii^a—O 

We will make the assumption that the expression in Eqn. [^may be reduced to the form 

OO 

Gi{z,t) = ^{z,t)(j){z,t) e{g{z,t))’‘^'^>’Pi-l{m^^b,T\mo) (10) 

nii^b—O 

where 9 and (f are functions to be determined, given knowledge of a particular partitioning rule. The partitioning mechanisms 
that we will subsequently consider can all be cast in this form, as we will demonstrate. 

We now consider the overall generating function describing a set of cell divisions. To represent the sum over all possible copy 
numbers before and after all cell divisions between divisions j and i, we introduce the notation 

/ OO P^iyb OO '^i — l,b OO P^j,b 

E- E E E E - E E 


( 11 ) 
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This combination of sums takes into account all possible states before and after each cell division i,i — for i > j. 

We note that the ordering of sums here progresses backwards in time from left to right: the leftmost two sums sum over all 
configurations related to the most recent cell division i, the next two sum over all configurations related to the preceding cell 
division i — 1, and so on. The final probability distribution is then 


/ n—1 

Pr,{m,t\mo) = '^Po{m,t\mn,a) 

n,l i—1 

where is a ‘probability propagator’ of the form 

= Ps{mi^a\mi^b)Poimi,b,T\m^_i^a), (13) 

representing the probability that a cell, which started with rrii-i^a units after division i — 1, grew to have rrii^b units, of which 
rrii a units were inherited by the next daughter cell after division i. The chain of divisions can be terminated at n divisions in 
the past by setting TOo,a = wq as the initial condition of the ancestor cell. Thus, a subscript 1 labels the first cell division, and 
a subscript n labels the most recent of n cell divisions. 

The overall generating function after n divisions is 


Gn{z,t\mo) 

/ 


i-1 


= y^y^2:""Tb(TO,t|mn,a) ]~[ 

m n,l 

/ n—1 


n,l i—1 

f OO T^n,b 

E E E ^ 

n— 1,1 mn^h — 0 ?7lTi,a=0 



(14) 


(15) 

n—2 


)n«* 

2=1 

(16) 


Generalising the approach of Rausenberger & Kollmann |31j . we now employ the assumption in Eqn. 10 to write 


/ OO Tnn,b 

E E E ^ 

n—1,1 rrin^h—O run^a—^ 

I 


n—2,1 mn,b—0 
I 

I E G(^(5(‘ 

n-2,1 

! 


Ps{mn,a\mn,b)Po{'rnn,b, T | m „_ 

n—2 

-l.a) 

(17) 

n-2 

Po{mn,b,T\mn-l,a) R 

2=1 

(18) 

2 = 1 

n—2 

\mn-l,a) 11 


(19) 

2=1 

n-2 

.) n 

2=1 


( 20 ) 


where in the Eqn. 20 we have changed variables zi = 6[g{z,t)). Comparing Eqns. 15 and 20 we can see that this process 
can be followed inductively. Each further step through the induction corresponds to another cell division, extracts a prefactor 
of (j){zi,T), and causes another change of variables Zi+i = 0{g(zi,T)). Extending this induction to n cell divisions, the overall 
generating function after n divisions is then 


G„{z,t) = ( 


\i=l 



^{z,t)ho{z,ty 


( 21 ) 


where hi and Zi are the solutions to 


hi{z,t) = g {e{h,+i),T) ; hn{z,t) = g{z,t) 
z^+i = e{g{zi,T)); zi = e{g{z,t)). 


( 22 ) 

(23) 


The recurrence relations Eqns. |22||23| are rather similar, but we retain their distinction for mathematical convenience, also 
noting that their indexing runs in opposite directions through time. Hence, the boundary condition for hi arises from the most 
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also arises from the most recent division, but in this 
21 vanish and the boundary condition hn = ho = g{z, t) leads to 


recent cell division, corresponding to i = n; the boundary condition for z, 
indexing this corresponds to * = 1 . 

It can be noted that for n = 0, the products in Eqn. 

Go{z,t) = as required. 

For BID dynamics, the product and final term in Eqn. [^are analytically tractable for several important inheritance regimes, 
allowing us to write the generating function in an exact form. We will first analyse two partitioning regimes of importance for 
biological modelling, and illustrate how this approach produces statistics of interest for mtDNA populations under these regimes. 
We will later explore other partitioning regimes. 


C. Binomial inheritance 

We first consider the case where agents are partitioned binomially at cell divisions. In this case, the following identities hold: 

Psimi^a\m^^b) = 

rrii^a—O 

and so (j){z, t) = 

where Eqns. [Mp7| follow by comparing Eqn. [^with Eqn. We are thus concerned with the solutions to the recurrence 
relations 




g{z,t)'^'’‘^ Ps{mi^a\mi^b) 




1 ; 


1 , 9{z,t) 
2^ 2 


(24) 


(25) 

(26) 
(27) 




hi 


1 , 9iZi-l,T) _ 1 

2 - 2 - 2 

( 1 I ^i+l 

n2+^'" 


hri. — 


9iz,t) 

2 


9{z,t). 


(28) 

(29) 


We will introduce the symbols I = and I' = for convenience here and throughout. In the Appendix we solve 

these related systems of equations, showing that the solutions take the form 


/Cii 2 ® + K 12 P 
nio2^ + Kuh ’ 
/v2i2® + K22^* 

^^ 232 * + ^24^* 


(30) 

(31) 


with Kii = l'^l'{z — 1)(A + 1/(1 — 2)), ki 2 = Ki 4 = 2"(A(/' — z{l + 1' — 2)) + iy{l — 2)), K 13 = Gl'{z — 1)A(^ — 1), K 21 = K 23 = 
1) + (Z — 2)(A2: — v)),K 22 = l'{z — 1)(?(A + 12 ) — 2v),K2i = 2X1'{z — 1){1 — 1). 

We also show in the Appendix that the first product in Eqn. takes the form 


Bn+iB-''-\A2 + B2){-A,/B,-pA/pB)u+i 

{Ai + Bi){ — A2/B 2 ', Pa!PB) n+l 

where Ai = 2X1'{I - l){z - 1), A 2 = XI'{z - l)(e(^-‘^)^(l - 2) + 1), Bi = B 2 = 1{XI'{1 - z) - {I - l)(Az - i/)), pA = hPB = 2, 
and 7 = a/A. (a; q)n is the g-Pochhammer symbol defined by (a; g)„ = ~ aq^). While this symbol is hard to interpret 

intuitively, we will see that expressions for important moments of distributions often only involve particular derivatives of the 
symbol that reduce to simple algebraic expressions (see Results). In addition, this term vanishes in the a = 0 case where 
immigration dynamics can be ignored. 

Combining ho from Eqn. [^and Eqn. we therefore use Eqn. [^to yield a closed-form generating function for BID dynamics 
with binomial partitioning at cell divisions (the full form is explicitly presented in the Appendix). In the Results section we will 
demonstrate the efficacy of this generating function solution (and subsequent solutions) by showing that moments derived from 
the generating function exactly match stochastic simulation (Fig. [^. This algebraic generating function for an arbitrary number 
of cell divisions extends a previous solution presented in integral form from Ref. |3I], with the advantage that the methodology 
allows straightforward analytic investigation of this and a wider class of systems; we next demonstrate this generalisation with 
a new inheritance regime which we will demonstrate has biological applicability. 















6 


D. Random additive or subtractive inheritance 

We now consider the case where a number of agents are lost or gained at each cell division, and this number is itself a random 
variable (we will consider the case where this number is a fixed constant later). For mathematical convenience we shall first 
assume that a certain number of agents are lost at partitioning, and that this number is taken from a binomial distribution with 
population size 2rj and probability so that the average loss number is rj (it is straightforward to see that a negative value 
for the 77 parameter corresponds to the gain of a number of elements identically distributed). In this case, by considering the 
possible values of n, the number of agents lost at a division, we obtain 


277 






n—0 


(33) 


0 

277 


= 0 + 


1 


1 


2 2g{z,t) 




and so (l){z,t) 


1 


1 


277 


,2 2g{z,t)^ 

9{z,t). 

The solutions to the corresponding recurrence relations are derived in the Appendix and are: 


(34) 

(35) 

(36) 


hi = 


Pl'iy{z-l) + h{iy-Xz)_ 
l^l'X{z-l) + h{ty- XzY 
ri'v{z — 1 ) + l{v — Xz) 
in'x{z-i) + i{v-xzy 


(37) 

(38) 


The first product in Eqn. 21 is ni=i?( Zi,T), which we show in the Appendix takes the form of Eqn. 32 with Ai = XI'{z — 


1), A 2 = XI'{z — 1), Bi = B 2 = l{iz—Xz), Pa = I, pb = l ,7 = a/A. The second product term is nr=i(l/2 + l/(2ff(-Zi, t)))^’^. 
We here introduce the symbols xi = X{l~^ — l)/{X — i'),X 2 = X{l — l)/{X — i'). In the Appendix we show that this product takes the 
form of Eqn. 32 with Ai = Z'(A+z/)(—X 2)"'(2 —1), = B 2 = — 2a:"(Az—i^), A 2 = 2l'v{—X2)'''{z — l),pA = xi,pB = (—X 2 ), 7 = 2g. 

We then have a closed-form generating function for random subtractive partitioning and BID dynamics (the full form is 
presented in the Appendix). We note here that applying this approach to the case of loss at cell divisions does not explicitly 
restrict copy number to be non-negative, and care must therefore be taken in its application. If the dynamics under investigation 
are such that the probability of copy number m < t] is negligible, the absence of this restriction will have negligible influence 
on results extracted from the analysis. If low copy numbers are likely, this approach can still yield useful results if a = 0, if 
expressions for P{m,t) are derived and P{m < 0,t) is treated as equivalent to P(m = 0,t)- This approximation is valid because 
the birth and death operations have m = 0 as an absorbing state, so a copy number below zero can never subsequently exceed 
m = 0. However, the simple expressions for E(m,t) and Y(rnA) in terms of generating function derivatives will yield incorrect 
results in these cases. The approximation will fail in cases where a non-negligible probability of attaining zero copy number is 
coupled with dynamics involving immigration (as opposed to birth). In this case, a specific boundary rule must be written in 
Eqn. [2 we have been unable to find closed-form solutions in this case. 


III. RESULTS AND APPLICATIONS 
A. Comparison with stochastic simulation 

In Fig. [^we compare the analytic results for copy number mean and variance, derived from our generating functions above, 
with the results obtained over 10® simulations using Gillespie’s stochastic simulation algorithm [32] • In order to compute 
moments arising from the generating functions for subtractive inheritance it is necessary to compute a small number of values 
corresponding to derivatives of the g-Pochhammer symbol (a; q)n with respect to the first argument a. We are not aware of a 
general analytic form for this derivative, but these values can be evaluated to arbitrary precision by symbolic software or through 
numerical perturbation, by evaluating ((a; g)„ — (a -I- e; q)n)/^ for suitably small e. In this case we take ‘suitably small’ to mean 
‘yielding a estimate converged to the required degree of accuracy’. In addition, in several cases which arise (for example, a = 0, 
emerging from our analysis below), this perturbative approach yields an analytic solution. In the Supplementary Information 
we provide Mathematica notebooks illustrating these calculations (and other calculations in this article). 

We choose arbitrary parameterisations for these confirmatory studies: we choose mo = 100, r = 10, for binomial partitioning 
we use a = 0.2, A = 0.06, iz = 0.01, and for subtractive inheritance we use a = 10, A = 0.01, v = 0.02, with 77 = 50 for deterministic 
and random subtraction. Analytic results exactly match simulation results throughout. 

We proceed by considering two examples motivated by specific biological questions involving cellular populations of mitochon¬ 
drial DNA (mtDNA). 
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FIG. 2: Comparison of analytic results for first two moments with stochastic simulation. Trajectories of copy number mean 
and standard deviation resulting from our analytic results (lines) and stochastic simulation (points). (A) Birth-immigration-death (BID) 
dynamics within cell cycles; binomial partitioning at cell divisions. (B) BID dynamics within cell cycles; loss of a binomially-distributed 
number of agents (p = = 100) at each cell division. Other parameters for these systems are given in the text. Theory and simulation 

match exactly in all cases. 


B. MtDNA bottlenecking: birth dynamics, binomial partitioning, changing population size 

In the case of the birth-death model with binomial partitioning, we can explore the levels of cellular noise introduced by 
controlled variation of the population size of cellular components without making continuous or steady-state approximations. 

We will consider a number rmax of dynamic phases labelled by r, where the rates A, v are constant within a phase but may 
take different values in different phases. To extend the above reasoning to describe different phases of dynamics it is necessary 
to compute the function gr{z,t) for each regime r, where gr{z,t) is the generating function using the appropriate parameters 
Xr, Vr for phase r and calculated at n = the number of cell cycles of phase r. For consistency with the above approach, we 
label phases starting from a zero index, so the first phase corresponds to r = 0, and we use Tmax to denote the label of the final 
phase. Then we use 



(39) 

hr = gr{hr-\-1^0^ 

(40) 

Goverall = ) 

(41) 


using induction over the different phases in the way way we used induction over different cell cycles above. Here we consider 
the changeover between phases by using the generating function at the start of the incoming phase. 

The system Eqns. 39|^ can be solved for arbitrarily many phases with different kinetic parameters, producing closed-form 
results for a wide range of different dynamic trajectories, including arbitrarily varying population size and cell doubling times 
(see Appendix). We illustrate this approach with the following simple two-phase model system. Initially, mo agents exist in a 
cell, with no associated cell-to-cell variability. These agents subsequently follow birth-only dynamics and binomial partitioning 
at cell divisions. The rate of birth is initially Ai in the first dynamic phase, changing to A 2 after ni cell divisions. We choose 
A2 = 21n2/T — Ai, to ensure that mean copy number returns to mo after a further ni cell divisions. We will use A2 > Ai < ln2/r, 
so the mean copy number either remains constant or initially drops to a minimum (the ‘bottleneck’) before recovering. 

For simplicity, we have here assumed that t , cell cycle length, is the same constant in each dynamic phase. This assumption 
can readily be relaxed by using different r^, so that cell cycle length is labelled by the current dynamic phase. In this way, our 
formalism can be used to explore the dynamics of systems with arbitrarily varying (though deterministic) cell cycle lengths. 

Using Eqn. 21 with Eqns. 30 and 32 the solution to Eqns. SOHS for two dynamic phases is simply G = gi{g 2 {z,t), 0)™“, with 
a = 0 and v = 0. After some algebra, we obtain, after ni cell divisions in the first phase and another ni in the second. 


V(m) 


moe-^i”i'"(e^i^ -k 2)(e^i"i^ - 2"^) 
e^i'^ - 2 


(42) 


The minimum copy number attained immediately follows cell division ni and is thus of size b = . This allows us 

to write the parameter Ai in terms of the bottleneck it produces, Ai = ln(2"i6/mo)/(?T,ir). Inserting this expression into Eqn. 
gives 
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FIG. 3: Modelling copy number variability due to mtDNA bottlenecking. (A) Trajectories of mean copy number of cellular 
agents born with rate Ai {t < 120) and A 2 = 21n2/r — Ai (t > 120). Lower Ai values enforce a smaller copy number bottleneck, with copy 
number recovering to its initial value after the bottleneck. Ai = In 2r is the value required to maintain constant mean copy number between 
cell divisions. Lines are analytic results; points are stochastic simulations. (B) Trajectories of coefficient of variation (CV) as different 
bottleneck sizes are imposed on the system. Horizontal fines give the analytic predictions for final CV derived from Eqn. |43| Other lines 
are analytic results; points are stochastic simulations. (C) CV as a function of bottleneck size from Eqn. 43 points show specific instances 
of stochastic simulation. 


_ mo(b- Too) / (6/mo)^/"i + 1 \ 

b \ — 1/ 


(43) 


Mitochondrial DNA (mtDNA) is observed to be present in fertilised oocytes at copy numbers around 10® [33l - l^ . During 
subsequent development, a pronounced copy number decrease occurs, as cells divide rapidly with little replication of mtDNA. 
The copy number per cell falls to a low bottleneck then recovers during later development. This mechanism is believed to 
ameliorate the inheritance of mutated mtDNA by increasing the cell-to-cell variability of mutant load in a cellular population 
and hence allowing cell-level selection to discard those cells that drift towards high mutant content. Substantial debate surrounds 
this topic [Si 137]: competing mechanisms have been proposed to increase mutant load variability [341135j , and the size of the 
copy number bottleneck, its power to generate variance, and thus its biological importance have been questioned [34] . 

Results from classical population genetics [Si |S9| have been applied to the statistics of mtDNA populations [40] , but even 
with refinements modelling fluctuations in the size of, and substructure in, the mtDNA population SDii, these results lack 
straightforward physical interpretability and the ability to address population statistics at arbitrary, non-steady-state points in 
developmental dynamics. The stochastic formalism we present has recently been used to address these issues, specifically in 
modelling the mtDNA bottleneck in mice [43] . Here we explore a more general question: what increase in copy number variance 
is possible through enforcing an mtDNA bottleneck of a specific size, and how does it relate to the size of that bottleneck? Results 
on the copy number statistics of mtDNA can then readily be extended to explore statistics of the mutant load of mtDNA, by 
consideriiw two decoupled populations of mutant and wildtype mtDNA. 

In Fig. [^we simulate the above system for ni = 12, mg = 10®, roughly matching the above mtDNA copy number magnitudes 
and number of divisions [44] observed in the mouse germ line. We use various Ai values, reporting mean copy number and copy 
number coefficient of variation (CV). It can be observed that lower Ai rates lead to more pronounced copy number bottlenecks, 
resulting in increased CVs that match the predictions derived from the above analysis. 

Eqn. [^ thus analytically describes the increased variance due to a copy number bottleneck under specific circumstances, 
whereby a copy number mg is decreased to a minimum b over ni cell divisions then raised to its original value over a further 
ni cell divisions. Lower bottleneck sizes b lead to exponential increases in the cell-to-cell variability associated with an mtDNA 
population. The assumptions within our illustrative model here can also straightforwardly be relaxed and closed-form solutions 
derived for more general dynamics, and a closed-form expression for the probability distribution function can also be derived 
using the above approach. We note than in this example, the variance does not converge to a final fixed value: longer times will 
result in higher variances. This feature can be altered using a model involving more homeostatic dynamics (see next subsection) 
or, in biology, may conceivably be dealt with on a cellular level by retaining cells with certain copy number statistics. 


C. Relaxed replication of mtDNA: immigration-death dynamics and various inheritance regimes 

A quantitative model for mtDNA dynamics throughout organismal lifetimes has been proposed to account for the intuitive 
feature that mtDNA copy number should be subject to cellular control [35] . This ‘relaxed replication’ model has influenced a 
wide range of studies on mtDNA dynamics in many contexts from human disease [46] to forensics [47] ; its quantitative behaviour 
has been explored (considering low-order moments without cell divisions) in the contexts of nuclear control on mtDNA [3S| and 
through simulation studies in topics including ageing mtDNA [49], the effect of anti-retroviral drugs on mitochondrial [50], and 
many others. As variability in mtDNA can have important physiological consequences [36] . we aim here to analytically extend 
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this model beyond a mean-only treatment to a more general (and realistic) situation both explicitly modelling the stochastic 
dynamics of individual mtDNA production and degradation and including different forms of cell division dynamics. 

The governing dynamics of copy number m in the model are 


dm amoptlr — ( 
dt \ —mlr 


i/t if to < 


— a —1 

otherwise. 


(44) 


where mopt is a ‘target’ copy number, f is the timescale of mtDNA degradation, and a > 1 is a parameter of the model describing 
nuclear feedback (in the original papers, a and f are respectively assigned the symbols a and r: we use tildes to avoid ambiguity 
with our analysis). For d not much greater than 1 and an initial condition toq < mopt, tbe probability of to > dmoptUct — 1) is 
very low; we will thus assume that the contribution of the term in Eqn. 44 corresponding to to > aTOopt/(a — 1) is negligible. 

When considering cell divisions in the above model, the meaning of mopt needs to be made explicit. We will consider mopt to 
be the target copy number for the end of a cell cycle, immediately before division. We then note that the first term in Eqn. 
contains a combination of a immigration term (independent of copy number) and a death term (dependent on copy number); we 
can therefore write the model as Eqn. [^with a = jdmopt, A = 0, and v = jd, defining the new parameter /3 = d/f. We can then 
use the above treatment to obtain the generating function of the relaxed replication model under different partitioning regimes 
(full forms given in the Appendix) and moments of interest of the copy number distribution in each case. 

Binomial partitioning. First we consider the case where mtDNA molecules are binomially partitioned at cell divisions. 
The dynamics in this case are illustrated in Fig. |^for mopt = m^ = 1000, d = 5, f = 10, and cell cycle length r = 5. This 
parameter set was chosen for compatibility with original work on the relaxed replication model: mtDNA copy numbers around 
1000 are biological reasonable, d = 5 is an intermediate value of the nuclear feedback parameter explored in Ref. [15], and 
f = 10 corresponds to an mtDNA half-life of 10 In 2 ~ 7, compatible with the range (in days) of half-lives assumed in Ref. [48] . 
The cell cycle length t = 5 was chosen both for rough biological applicability (corresponding to cells dividing every 5 days) and 
to illustrate transient behaviour of the model. 

It can be observed that the variance and the mean converge on the same behaviour; the difference between the two is 
straightforwardly found to be E(to) — V(to) = 4“"TOoe“^^*^*“''”'’’\ clearly decreasing to zero with cell divisions n. The mean 
copy number immediately before a division E(to,t) takes a value that approaches, but does not reach, mopt- Some algebra (see 
Appendix) shows that this value, in the long time limit, is 


E(TO,r) = mopt - 2g/3T _ ^ '^opt- 


(45) 


The generating function analysis also allows an expression to be derived for the probability distribution function for mtDNA 
copy number (see Appendix). 


P{m,t) = (l/TO!)(-6)’"(l-&)’”'>-'"e-“C/(-TO,l-TO-f TOo,a(5- l)/6), (46) 

where [/(.,.,.) is the confluent hypergeometric function, a = TOopt(l — -f (2“”e“^*(e^'’’ — 1)(2” — e“^"'’’))/(2e^'’’ — 1)) and 
b = e ~As n —>■ oo, this converges to the simpler expression 


P{m,t) = 


^ TOopt(l-2e^" + e-^(*-")) ^^^ TOopt(l-2e^"+e-^(*-")) 


m 




2ed'^ - 1 


2eP'^ - 1 


(47) 


Eqn. 46 thus provides a complete solution for the relaxed replication model with binomial cell divisions (see Fig. for a 
comparison with stochastic simulation for E(to), V(to) and P{m); other statistics and quantities of interest are readily extracted 
from the generating function). 

Subtractive partitioning. Next we consider the case of dynamics under which a random amount of mtDNA is lost at 
each division. This picture could model, for example, mtDNA dynamics in budding yeast, where cells with 50 — 200 mtDNA 
molecules m undergo asymmetric partitioning, with 10 — 20% of their mitochondrial content being lost at budding events, and 
homeostasis acting to maintain copy number |52j . Fig. |4 illustrates the behaviour of this system with mopt = 100, r; = 15. 
Interestingly in this case, the variance of copy number reaches a minimum at an intermediate point in each cell cycle, after the 
partitioning event (which increases the variance) and before an extended period of dynamics under homeostasis has led to an 
increase in copy number and variance. As above, the difference between mopt and E(to,t) can easily be computed in the long 
time limit (see Appendix): 


E(to, r) = mopt - ^ ■ (48) 

The expression for the variance in this case involves derivatives of the q-Pochhammer symbol, which at first would seem to 
prevent extracting simple and intuitive expressions for the variance. However, as with several applications of this approach, the 
derivatives involved all reduce to simple algebraic expressions (see Appendix). Taking the limit of many cell divisions n ^ oo, 
we obtain 
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FIG. 4: Copy number statistics for the relaxed replication model of mtDNA. Lines are analytic results; points are from stochastic simulation. 
(A) E(m) and (B) V(m) for binomial partitioning, modelling dividing cells, with mo = niopt = 1000, /? = |. V(m) and E(m) converge to 
the same trajectory; grey line gives analytic result for E(m) at the end of a cell cycle. (C) E(m) and (D) V(m) for random subtractive 
partitioning, modelling budding cells, with mo = 100 = rriopt = 100, /3 = |, r; = 15. Grey lines give analytic results for E(m) at the 
end of a cell cycle, and V(m) at start, end, and minimum point (blue crosses) in the cell cycle. (E) Comparison of theory and stochastic 
simulation in illustrative snapshots of probability distribution functions at given times in the binomial case using Eqn. |46| 


V(m, t) = (3e2/5(-‘) - . (49) 

This expression allows us to characterise the form of the variance curve. Writing r]' = -q/ (2(e^^’’ — 1)), we find that V(m, 0) = 
mopt + Y{m,T) = mopt + q'{^ — and that the minimum variance occurs at f' = t — l//31n (|(1 + 

and takes value V(m,t') = mopt — {q/6) coth.{f3T/2) . All these results agree exactly with stochastic simulation, as illustrated in 
Fig.H and other statistics and quantities of interest can readily be extracted from the generating function. 

We note that as this example falls in the regime where q <C mopt, so the probability of low copy number is negligible, 
the statistics derived using our generating function approach are reliable (and the approach can also be used to compute the 
probability distribution function and other moments). In cases where low copy numbers are likely, caution must be taken in 
employing this approach, as described above. 


Extinction probabilities under balanced copy number dynamics 

A strength of the use of generating functions to analyse stochastic dynamics and partitioning of cellular species is that statistics 
other than low-order moments can be straightforwardly computed. As demonstrated in the previous subsection, full probability 
distribution functions can be extracted for cellular populations of agents from generating functions, although these functions 
can be rather complicated. As a simpler example of biological interest, we here consider the extinction probability P{0,t) of a 
cellular species, under the two dynamic regimes we have previously considered in the context of mtDNA dynamics. The first 
example is birth-death dynamics with binomial partitioning, with A = In 2/t -|- as used in the bottleneck section. The second 
is immigration-death dynamics with binomial partitioning, as used in the relaxed replication section. Both of these examples 
exhibit a balanced mean copy number, with the expected production of agents over a cell cycle balancing the expected loss 
through cell divisions. As previously described, the variance of the birth-death case increases with time, whereas the variance of 
the immigration-death case converges. 

The fact that extinction probability can be straightforwardly extracted from onr generating functions, using P{m = 0) = 
G'(0,t), allows us to explore the probability of extinction under these dynamics. The resulting expressions under birth-death 
(BD) and immigration-only (I) models are 


Pbd{ 0, t) 


Pi{0,t) 


2*/'^(^T(n -f 2) -I- n In 2) — 2vt 
2*/"^ (z/r(n -f 2) + (n -I- 2) In 2) — 2z/r ^ 

mopt (2 -"(m - l)u" - (u - 2)u' 


(1 - 2-”u"u')™“ exp 


u'{u — 2) 



(50) 

(51) 


where u = e~^'^ and u' = Consideration of the n —>■ oo limit shows that in the long time limit, extinction probability under 

the birth-death model converges to unity, whereas in the immigration-only model a limiting probability is reached. Setting t = 0 
for simplicity (thus considering the population at the start of a cell cycle), this limiting probability is exp(—mopt(l -f (u — 2)“^)). 
The difference between the ID and B cases is due to the irreversibility of extinction under the birth-death model (and the 
non-zero probability associated with extinction during every cell cycle); by contrast, extinction in the immigration model can 
be escaped as immigration creates more agents in the cell withont necessitating a nonzero source population. Hence, a nonzero 
probability flux away from the m = 0 state exists, and eventually balances the flux into that state due to partitioning noise: 
the extinction probability may thus be thought of as representing the proportion of time during which the system occupies the 
m = 0 state. 

Further results can straightforwardly be extracted from our formalism for extinction probabilities in non-balanced cases. As 
a brief example, we consider the birth-death model with binomial partitioning, with a new parameter k = \ — v — In2/T, so 
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FIG. 5: Solutions from recurrence relations for different partitioning regimes. Points are from stochastic simulation, lines 
are theoretical predictions. Deterministic, clustered, subtractive and binomial inheritance regimes are compared. Parameters used were 
mo = 200,1^ = 0.01, T = 10, ric = 10, A = log2/T + u,ri = mo - the choice of the latter two parameters was fixed to preserve a constant 
mean copy number. Stochastic simulation results are from an ensemble of 10® simulations of each situation. 


that K measures the ‘excess’ birth rate beyond that required for copy number balance. For k < 0, extinction is certain in the 
long-time limit, but for k > 0 there is a finite probability that the copy number will never reach zero. To illustrate the qualitative 
behaviour of the system, we set = 0, and we obtain 


PBD{m = 0 ) 


1 - e'' 


\ _|_ ^KTIT _ 


(52) 


showing that extinction probability in the long-time limit decays roughly exponentially with k. The more general extinction 
probabilities for u ^ 0 and in the absence of cell divisions are given in the Appendix. 


Other inheritance dynamics 

In addition to the binomial partitioning and random additive or subtractive changes of copy number, we have explored several 
other possible dynamic regimes of inheritance. The case where a fixed, deterministic number of agents is gained or lost at cell 
divisions is analytically tractable (see Appendix). We also consider deterministic halving of copy number, so that each daughter 
cell inherits half of the mother cell’s content (rounded down). Additionally, we consider the inheritance of clusters of agents, such 
that agents are split into clusters of size ric and these clusters are binomially partitioned. We have not found closed-form analytic 
solutions for the generating function for a general number of cell divisions for these latter two cases, but analytic statistics can 
nonetheless be obtained for a given number of cell divisions through the calculation of the appropriate recurrence relations. 

Fig. 0 illustrates the use of this approach to calculate the mean and variance of copy number for these systems, and for 
birth-death dynamics in the binomial and constant subtractive inheritance regimes described earlier. The agreement between 
stochastic simulation and analytic results is again excellent, showing, as expected, that deterministic inheritance leads to the 
lowest magnitude of stochasticity in copy number, followed by binomial partitioning, followed by clustered partitioning (illustrated 
for Tic = 10 in this case). We expect that other inheritance regimes of interest may be addressable through a similar approach. 


IV. DISCUSSION 

We have introduced a general mathematical formalism with which to address the stochastic dynamics of cellular agents that 
are inherited according to non-trivial and potentially stochastic dynamics at cell divisions. Our approach differs from, and 
extends, several previous tools designed to address stochastic partitioning in biology. First, our approach yields full, closed-form 
generating functions for several cases, allowing the extraction of all details of copy number distributions, rather than focussing 
on variance or other low-order moments alone. Second, we nowhere assume that a steady-state or equilibrium has been reached, 
and are thus capable of extracting copy number statistics at any given time during the stochastic biological process of interest. 
Third, we focus on birth-immigration-death dynamics rather than the more common stochastic gene expression dynamics, with a 
view to modelling the behaviour of non-protein cellular components (including mtDNA, which we explore in particular). Fourth, 
we explore several specific partitioning regimes, obtaining closed-form results for arbitrary numbers of cell divisions under those 
we term binomial partitioning and random and deterministic subtractive or additive inheritance. We also obtain results for finite 
numbers of cell divisions under deterministic partitioning and binomial partitioning of clusters. 

We have focussed on agents undergoing birth-immigration-death (BID) dynamics between cell divisions; we expect that any 
random dynamics with a corresponding generating function of the form Eqn. will also admit treatment using this approach, 
paving the way for further generalisation of this approach. 
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We note that the results arising from our analysis should be interpreted as ensemble statistics of single-cell measurements. If 
tracking a particular lineage of dividing cells, it should be remembered that the statistics of daughter cells will exhibit correlations 
due to their shared heritage. Our results represent expected statistics from a well-mixed bulk case. 

In the absence of immigration dynamics, we have shown that a system governed by different rate parameters at different times 
can also be solved analytically. The rates associated with the production and destruction of cellular agents can vary arbitrarily 
as long as the rate of this change is lower than the cell division rate. This approach therefore provides a way to explore the 
statistics of stochastic systems with arbitrarily-changing population size, by contrast with many results from classical statistical 
genetics which assume a constant or constant-mean population size (and additionally often only provide equilbrium results, 
preventing the quantitative exploration of systems before a steady state is reached) [351 IMl SIl ■ It is also straightforward to 
vary the cell cycle length r in these dynamic phases and so allow a treatment of cellular dynamics under varying division times. 
We have illustrated a general use of this approach in addressing the mtDNA bottleneck, and it has been used to explore the 
bottleneck specifically in mice in detail [43] . We believe that this formalism may prove useful in other contexts where organellar 
content is subject to dramatic and non-random population size changes, for example, in considering cellular populations during 
tumour development, where variability in cellular conditions causes time differences in physiological rate constants as tumour 
cells continually divide [53j . 

We have demonstrated the applicability of our stochastic formalism with some illustrative problems from cellular biology. We 
have found expressions for the cell-to-cell variability in mtDNA populations due to the imposition of a copy number bottleneck 
of given size, and extended a the classic ‘relaxed replication’ model of mtDNA to include the stochastic dynamics of individual 
mtDNAs, and the effects of cell divisions. This model is widely influential in the study of mtDNA genetics and disease, but its 
quantitative analysis has typically been limited to descriptions of its mean behaviour or simulation studies focussed on used. 
We have thus used our approach to further analytic understanding of this important model. Furthermore, we have explored in 
detail the statistics of populations of cellular agents under passive copy number balance over many cell divisions, a situation of 
importance for organelles and which may be of general applicability in cell biology. 

Accurate models for the variability of cellular populations enable more powerful inference using experimental measurements of 
mean and variance across cells |S3] . We hope that our results for stochastic inheritance dynamics will facilitate the strengthening 
of this link between theoretical and experimental biology and allow more information about underlying cellular dynamics to be 
obtained from the wealth of experimental measurements currently appearing. 
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Appendix 

In this Appendix we present some of the lengthier mathematical results used in the main text. We include a Mathematica 
notebook as part of this article, containing the following derivations. 

Generating function for the BID process 

The derivation of the generating function for BID dynamics in the absence of cell divisions is well known: we include it here 
for completeness. We write the PDE describing the generating function G{z, t) in Laplace form: 


dG{z, t) 
dt 

- {iy{l - z) + X{z'^ - = a{z - 1)G 

(53) 

G(z,0) = z^^, 

(54) 



(55) 


We proceed by using the method of characteristics, writing down ODEs describing how the parameters of G, and G itself, 
changes along a characteristic curve, with progress along such a curve parameterised by s. The corresponding ODEs are 


dt 

ds 

dz 

ds 

dG 

ds 


1 


-{v{l - z) + A(z^ - z)) 
a{z — 1)G 


(56) 

(57) 

(58) 


Eqn. 56 lets us immediately set t = s, omitting a constant of integration as the absolute value of progress along a characteristic 
curve is unimportant. Using t = s throughout, Eqn. I^is solved by 


1 _ 

1 _ _\gCi(A-y)-t(A-!y) 


(59) 


where ci is a constant of integration, the explicit form of which will be useful later. Rearranging this into an expression for ci 
gives 


Einally, Eqn. 58 with Eqn. 59 gives us 


Cl =t + 



\ — V 


(60) 


G = c2e-“‘ . (61) 

C 2 is a function of ci because the quantity ci, the constant of integration acquired when integrating z with respect to 
s, is independent of s, and hence forms an independent parameter when integrating G with respect to s. We require that 
G{t = 0) = 2 ;"*“, so we choose 


C2(ci) = - Ae^^^) 


;,^-a/A - I 


VAeA-‘^)ci _ 1 


(62) 


where the first term cancels the final term in Eqn. [^when t = 0, and the final term can be seen to extract a factor z'^° from 
Eqn. I^for ci when t = 0. We then have 


G{z,t) = C 2 (ci(z,t))e-“* 

which, after inserting Eqn. [^and some algebra, gives 

/ v-\ f ve'^^-''">\z - 1) - Xz + v \ 

\ — 1) — Xz -\- vj \ Ae('^“‘')*(z — 1) — Xz + i') 


G{z,t) 


(64) 

(65) 
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Recurrence relations arising from induction over cell divisions 

This section focusses on the solution of recurrence relations of the form 

aCi+i +b ■ 1 ., A dQ-b 

Q = — -or equivalently, Q+i = —— (66) 

cQ+i + a -cQi + a 

In the Main Text, both hi and Zi follow relationships of this kind; we use the symbol C,i here to emphasise that the same 
solution strategy applies in both cases, and describe specific solutions below. This system is solved, after [55], by defining 
a = , P = ^, D = ad — be, yi = (i + d and implicitly defining Wi through yi = . These changes of variables allow us to 

find an expression for Wi, which can then be substituted back through the above chain to find Q. We have 

fd 

Vl = OL - 

y^+l 

/3Wj+2 

- = a - 

W^+l Wi+i 

-)> l3wi+2 - aWi+i + Wj = 0, 

which is solved by considering solutions to the characteristic equation j3k‘^ — afc + 1 
± — 4/1). Then 


(67) 

( 68 ) 
(69) 

= 0 , which are straightforwardly ki ,2 = 


Wi 

Vi 


0 


Cik\+C2k^ 

Cok\ + /c2 
Cokl+^ + kl+^ 

Cok\ + ^2 d 
Cofci+I + kl+^ ~ C 


(70) 

(71) 

(72) 


where Ci are constants to be determined from boundary conditions. If the boundary condition takes the form = yjpf, as 
is the case throughout the situations we consider, we obtain 


Cok\ k2 d 
Cofci+I + kl+^ ~ C 

^Co 


pz + q 
rz + s 

^2 fcr”(^ 2 c(pz + 9 ) + k 2 d{rz + s) — c{rz + s)) 
c(rz + s) — kic(pz + q) — kid(rz + s) 


(73) 

(74) 


Thus, given knowledge of a,b,c,d from the recurrence relation and p,q,r,s from the initial condition, we can obtain ki,k 2 
through a and /3 and hence use Eqns. 74 and 72 to obtain a solution to the recurrence relation. Below, we use this approach to 
obtain solutions for the systems of interest in the main text. 


Binomial partitioning solution for h 

We will use the substitutions I = I' = The original recurrence relation is 



{vl — X)hi+i + (—A — i/(Z — 2)) 
{X(l — l))/ii_|_i + (2^ — X{1 + 1)) 
hn = g{z,t) 

{vV — X)z + ( 1 / — vV) 

{X{1' - l))z + {v - Xl'Y 


(75) 

(76) 

(77) 

(78) 


hence, a = {ul — A), b = {2iy — vl — X),c= {XI — X),d= {2iy — X — Xl),p= {vl' — X),q= {v — vl'),r = {X{V — 1)), s = {v — XI'). 
Using these values to determine a,ki,k 2 , and after some algebra, we obtain 


2*rr(z - 1)(A + y{l - 2)) + 2"f (A(Z' - z{l + 1'- 2)) + v{l - 2)) 
2HH'X{1 - l){z - 1) + 2^h{X{l' - z{l + 1'- 2)) + v{l - 2)) 


( 79 ) 
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Subtractive partitioning solution for h 


g{hi+i,T) 

(80) 

(vl — X)hi+i + (u — vl) 

{X{l-l))h,+i + {iy-Xiy 

(81) 

9{z,t), 

(82) 

{vV — X)z + {v — vV) 

{x{i'- i))z + {v - xvy 

(83) 


hence a = (vl — A), b = {v — vl), c = {\{l — 1)), d = {v — XI),p = {uV — X),q = {v — vV), r = X{1' — 1), s = {v — XI'). Then 


l'^l'v{z-l) + l^{v-Xz) 
l^l'X{z — 1) + l'{v — Xz) 


Binomial partitioning solution for a 


(84) 


^ 1 g(z^-l,T) 

2~^ 2 

{1{X + ly) — 2X)zi-i + {2i> — 1{X + v)) 
{2X{l-l))z^-l + {2v-2Xl) 

(2z/ — 2X1)Zi-^i (^(A v) — 2v) 
(2A(1 — Z))zi+i + (^A + v)) 

1 , 9{z,t) 

2 ^ 2 ■ 

{l'{X + v)-2X)z + {2v-l'{X + v)) 
{2X{1' -l))z + {2v-2XV) ■ 


(85) 

( 86 ) 

(87) 

( 88 ) 

(89) 


In Eqn. [^we have used the equivalence in Eqn. [^to rewrite the recurrence relation in the form we have previously solved. 
Subsequently, a = {2v — 2X1), b = {1{X + ly) — 2v), c = (2A(1 — 1)), d = (Z(A + iy)),p = l'{X + v) — 2X),q = {2v — l'{X + iy)),r = 
{2X{1' — 1)), s = {2iy — 2X1'), leading to 


ri'{z - 1)(/(A + iy)- 2u) - 2 H{XI'{z -!) + (/- 2)(Az - v)) 
2XlH'{z - 1){1 - 1) - 2H{Xl'{z -1) + {1- 2){Xz - ly)) 


(90) 


Subtractive partitioning solution for 2 


Zi = g(z,_i,r), 

{lyl — X)zi-i + {ly — vl) 

{X(l — l))zi-i + {ly — XI) 
_ {ly — Xl)zi+i + [vl — ly) 

~ {X{l-l))z + {iyl-X) 

zi = g{z,t), 

[lyV — A)z + {v — lyl') 
{X{1' - l))z + {v - XI')’ 


(91) 

(92) 

(93) 

(94) 

(95) 


Hence a = {ly—Xl), b = {vl — v),c= {X{l — l)),d= {vl — X),p = V{X+v) — 2X), q = {2iy—l'{X+iy)), r = {2X{1' — 1)), s = {2v—2Xl'), 
leading to 


ri'v{z — 1) + — Xz) 

in'X{z-l) + l{iy- Xz) 


( 96 ) 
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Products of prefactors over the inductive process 


We are concerned with an expression for the product in Eqn. 21 ni=i We first consider ni=i ''') which occurs 


as a factor in this expression in every partitioning regime. We recall that has the form 




-A 


yg(A _ X) _ _\2 _|_ i 


Ot/\ 


(97) 


It will be convenient to write 


Zi 


S,{Z^,T) 


^iPa + ^iPb 

Mp\ + 

/ Aip\ + 

\A2P\ + B2p'‘g) 


(98) 

(99) 


where Ai — (^A 2 (y — A)), Bi — (^B 2 {v — A)), A 2 — (A/(Ai — A 2 ) + VA 2 — AAi), B 2 — (AZ(i?i — B 2 ) + VB 2 — Ai?i), 7 — <a/A, 
following from the form of ^(z,t). 

The product of an expression of this form can be written 




i—1 2=1 ^ ^ ' \ 2=1 

Here we make use of the g-Pochhammer symbol (a; g)„, dehned by 


1 + ^^ 
Bi Ph 




2=1 


{a-,q)n = “ “9^)- 


A;=0 


( 100 ) 


( 101 ) 


It can then be seen, by setting a = —Aj(Bj and q = pa/Pb, 


By’sU 

2=1 



1 + {-A,/Bf,pA/pB)n+l, 

Tjn+1 Pa/PB) n+l 

Pb -A. , R.- 


( 102 ) 

(103) 


and hence that 




2=1 


Mp^ 

Bi Pb 




2=1 


AI 2 Pa 
B 2 Pb 


B/+^B^^-\A2+B2)i-A,/B^;pA/pB)n+l 

(Hi + i3i)(—H 2 / B 2 ] Pa/PB) n+i 


(104) 


yielding a simple form for the product of interest. For the binomial case, we have from Eqn. 90 


-Zi = 


lH'{z - 1)(;(A + u)- 2v) - 2 H{XI'{z -1) + {1- 2){Xz - v)) 
2XlH'{z - 1){1 - 1) - 2%Xl'iz - 1) + (1 - 2)(Az - v)) 


(105) 


so that Hi = (l'(z-l)(Z(A + z/)-2z/)),Hi = {-l{Xl'{z-l) + {l- 2 ){Xz-,y))),A 2 = (2Al'(z-1)(1-1)),H 2 = (-;(Al'(z-1) + (1- 
2){Xz-i^))), hence Hi = 2XP{1 - l){z - l)(j 2 - A), H 2 = XW{1 - l)(z - l)(i^- A), Bi = B 2 = -l{XV{z -!) + (/- 2){Xz - v)){v-X), 
and PA = I, Pb = ‘2,"f = a/X. 


For the subtractive case, we have from Eqn. 96 


ri'v{z — 1) + l{y — Xz) 
lH'X{z-l) + l{v-Xz) 


(106) 


so that Hi = {l'v{z — \)),Bi = {l{v — Az)),H 2 = {l'X{z — 1)),H2 = ~ ^z)), hence Hi = Xl'{z — l){v — A),H 2 = 

Xll'{z — l)(i/ — A), Bi = B 2 = l{v — A)(i/ — Xz), and pA = I, Pb = = oi/X. 
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Other products 


We can also use this result to compute the product of exponentiated prefactors involved in the subtractive inheritance regimes. 
The first product required is HlLi We recall the definitions of the recurrence relations in the Main Text 


hi{z,t) = g{ 9 {hi+i),T) ] hn{z,t) = g{z,t) 

Zt+i = 0{9{zi,T)); zi = e{g{z,t)). 

For both subtractive inheritance cases, 9 {g{z,t)) = g{z,t), so it can straightforwardly be seen that 


(107) 

(108) 


Zi+i = g{zi,T ); zi = g{z,t) 
hi{z,t) = g{hi+i,T)-, hn{z,t) = g{z,t) 
and so/i„(z,t) = Zi] = Zi- 


(109) 

( 110 ) 
( 111 ) 


Thus, the product Jir^i '’’) ^ i® equivalent to the product nj=i where j = n — i + 1 . As * and j are dummy variables, 
we can then identify the required solution as rir=i We have from Eqn. 


84 


that 


riV(z-l)+Z*(z/-Az) 

+ Xz) 

iyl'x\{-X2T{z - 1) + x^{-X 2 y{iy - Xz) 
Xl'x\{ — X2)^{z — 1) + Xi{ — X2Yiv — Xz) 


( 112 ) 

(113) 


where we have rewritten the final line to 


avoid a diverging factor of (z — 1) ^ 


appearing in the Pochhammer symbol, using 


= ^(e(-^k-l) 
A — u 

X — u 


(114) 

(115) 


We then see that ^ is of the form Eqn. 99 so we can use the result therein, with Ai = vV{—X2)'^{z — l),i?i = B2 = 


x^{v - Xz), A2 = Xl'(-X2)"'(z - 1), pA = Xi,pB = (- 2 ^ 2 ), 7 = -V- 


Finally, we wish to compute an expression for Y\Y=i ( ^ 2 g{z t) ) ’ random subtractive regime. We again exploit 

to show that the kernel of the desired product is equivalent to Again, 


the relation between g{zi, t) and hj in Eqn. 
we use hi from Eqn. 1^ after some algebra 


this expression reduces to 


Vx\{—X2)^{z — 1)(A — v) -\- 2 xi{—X 2 y{v — Xz) 

21'x\{—X2)'^i'{z — 1 ) + 2 a :"(— X2 )®(^ ~ 


(116) 


whereupon we can use Eqn. 
xi,pb = (-2:2),7 = 277. 


104 with Ai 


V{-X2)'^{z - 1)(A - v), Bi = B2 = 2a;"(j/ - Az), A 2 = 21 '{-X2)'^v{z - 1), pA = 


Full forms of generating functions 

To recap, we use a, A, v to respectively represent the rates of immigration, birth, and death in our model; uiq for initial copy 
number; r for cell cycle length, n for the number of divisions that have occurred, and t for the elapsed time since the most recent 
cell division. We employ simplifying symbols I = V = and xi = X{ 1 ~^ — 1)/(A — t^); 2:2 = X {1 — 1)/(A — 12 ). 

The general form of the generating functions was shown in the Main Text to be 



(i) (ii) (iii) (iv) 


(117) 


where Zi and hi are the solutions to recursion relations defined in the Main Text. The term (iii) is the same in all calculations 
and is, from Eqn. |64[ 
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u — X 


Xl(^z — 1) — Xz i 


y./\ 


(118) 


The generating function in the case of binomial partitioning at cell divisions involves (i) Eqn. 104 applied to the appropriate 
terms described in Eqn. 105 (ii) unity; (iii) Eqn. 118 and (iv) from Eqn. 79 giving overall 


Gn{z,t) = 


' {Xll'{l - l){z - 1) - l{Xl'{z -1) + {1- 2){Xz - ly))) ( 


-2Xl'{l-l){z-l) . I 


y/X 


n+1 


^ {2X1'{I - l)(z - 1) - l{Xl'{z -1) + {1- 2){Xz - :/))) ( 


-Ai/'(/-l)(z-l) . I 

-i(A;'(z-l) + (i-2)(A2-i/)) > 2 


n+1 y 


-A 


a/X 


Xl{z — 1) — Xz + V ^ 

l''l'{z - 1)(A + y{l - 2)) + 2”(A(r - z{l + I' - 2)) + u{l - 2)) ^ 
I'^l'iz - 1)X{1 - 1) + 2^{X{1' - z{l + 1' - 2)) + v{l - 2)) 


(119) 


The generating function in the case of random subtractive inheritance involves (i) Eqn. 104 applied to the appropriate terms 
described in Eqn. 106[ (ii) Eqn. 104 applied to Eqn. 116 as described; (iii) Eqn. 118 and (iv) from Eqn. 84 and is 


Gn{z,t) = 


{21'v{-X2T{z - 1) - 2x?(Az - v)) 


)(-3:2)"(z-l) . XI 


— {Xz — iy) 


27? 


71 + 1 


{l'{X + v){-X 2 )''{z - 1) - 2x+Az - v)) 

\xi'l{z-l) + l{r,-Xz))(=§^-l) 


-X2Y{z-1) . XI 
— 2x'^{Xz — y) ’ —X 2 


71+1 , 


Ol/X 


71+1 


^XV{z-l) + l{,-Xz))(=^;^-l)^^^ 

v-X (l''l'v{z-l) + v-Xz 


^Xl{z — 1) — Xz + v 


l'^l'X{z — 1) + v — Xz 


( 120 ) 


Birth-death dynamics for the mtDNA bottleneck 


For birth-death dynamics, a = 0, so for binomial partitioning the generating function takes the form of the final term in Eqn. 

[m 

In the case of balanced copy number, with A = In2/T -|- z/. 


(2‘/^ - l)(z- 1) — z In 2 

^ ’ 2*/'^(z — l)(z/r-|-ln2) — zln2 — z 2 t(z — 1) 

The corresponding solutions for Zi, hi are 

2*/'^(z — l)((z -I- l)vT -f i In2) — 2(z — V)vt — z ln4 

* 2*/'^(z — l)(i -|- l)(zAT -|- ln2) — 2(z — l)z^r — zln4 

^ 2*/'^(z — l)((z — n — 2)ut + {i — n) ln2) -|- 2i>t{z — 1) -|- zln4 

* 2*/'^(z — l)(z — n — 2){vt -|- In2) -|- 2vt{z — 1) -|- z ln4 
so 


Gn{z,t) 


1 X 
(i) 


1 X 1 X 



2vt{z — 1) — 2‘/'^(z — l)((n -I- 2)vt + nln2) -f zln4 
2vt{z — 1) — 2*/'^(z — l)(n -l- 2,){vt -f ln2) -|- zln4 

■v' 

{iv) 


rriQ 


( 121 ) 


( 122 ) 

(123) 


(124) 


Relaxed replication of mtDNA dynamics 

The generating function for a single cell cycle (involving immigration and death dynamics) can straightforwardly be found 


G{z,t) = ® mopt{z - 1) {l + e ^‘(z-l))"*”. 


(125) 
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and so 


® nioptiz - 1) (126) 

g{z,t) = l + e-^\z-l). (127) 

For binomial partitioning, we have (j){z,t) = 1 and 9{g{z,t)) = (1/2 + g{z,t)/2). The solutions of recurrence relations Eqns. 
[22]IMl are then 


h, = l + 2*-”e-^(‘+("-*)^)(z-l) 


So the overall generating function is 


Gn{z,t) = exp ( mopt{z - 1) (1 - e + 


2-„g-/3t(g/3T _ _ e-/5"^) 

26/5-^ - 1 


(i) 


X (l + 2-”e-^(‘+”^)(z - 1))™“, 


which can be written as 


with 


(iv) 


G(z,t)=e‘^^+\cz + d)^°, 


(128) 

(129) 


xlxexp(l-e ^*mopt{z-l)) 


(ii) 


(130) 


(131) 


d — TTlopi 


1 - e-^‘ + 


2-„g-/3t(g/3T _ - e"^”"") 

2e/5^ - 1 


b = —a; 
d = 1 — c. 

In particular, the mean copy number is (ae“^+^(c 2 + d)™° + e“^+*’moc(cz + giving 


(132) 

(133) 

(134) 

(135) 


E(to, t) = 2-”e-'^(*+"^)mo + m^p* ( 1 - + 

and setting t = r (at the end of a cell cycle), we obtain 


2-ng-/3t(g/3r _ 2)(2" _ 

' 2e/5^ - 1 


(136) 


1 -I-^ , 

E(m,r) - ruopt = - ^ - -mopt - 2“”e“^'^(”+^)TOo. 


The n —> oo limit of this expression is 


, 1 -,/ \ n—foo 1 

E(m,T) - rUopt - 2e/3r _ 


In this n —>■ oo limit the generating function reduces to 


G{z, t) exp {mopt{z - 1) ^1 - e + 


g-/3t(-g-/3t _ 

2ef^'^ - 1 


so that 


(137) 


(138) 


(139) 


1 /mopt(l-2e^^ + e-'5(‘-^))\™ _ 2e^^ + 


m\ V — 1 


exp 




I 


(140) 
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Using Eqn. |131|and Leibniz’s rule, the general probability distribution function is given by 


X 1 

P{m,t) = — j— 
m\ oz'^ 

1 


z=0 


m 


d 


^m-k^az+b ™0- 

ml ^ \ k J (too — k)\ 


(cz + 


= —y 

to! ^ 

k=0 

e“^+'’(-c)™(cz + d)“« 

to ! 

= (1/to!)(—(—to, 1 — to + toq, —ad/c) 


U —TO, 1 — TO + Too, 


z=0 

—a{cz + d) 


2=0 


(141) 

(142) 

(143) 

(144) 


where U (a, 6, z) is the confluent hypergeometric function. 

For subtractive partitioning, we have (j){z,t) = (1/2 + l/(2(7(z, t)))^’' and 6{g{z,t)) = g{z,t). The solutions of recurrence 
relations Eqns. 22p3 are then 


So the overall generating function is 


z, = l + e-^(*+(*-i)")(z- 1 ) 
hi = l + e-^(*+("-*)^)(z- 1 ) 


(145) 

(146) 


G„(z,t) = exp(TO„pt(z-l)(e-''‘-e-'5(‘+"")))4’' 


' (z - 1 + e^(‘+”^)) 

(^ — 1 + 2e^(*+”'’‘l) e“d(*+"T)(z-i); 


2r, 


(i) 


n+1 


n+1 , 


(iz) 


( \ ^ 
l + e-^(‘+”")(z-l)j 


(147) 


(zu) 


(zu) 


In particular. 


E(to, t) = TTiopt + e (too - niopt + ??(1 + (0; 

It follows straightforwardly from the definition of the q-Pochhammer symbol that 


(148) 


(0|9)zz+i — 


qn+l _ I 


(0;q)"+i = q 


q-1 

{q^ - l)(g”+i - I) 


(g + l)(q-l)2 

Using these results and setting t = t (at the end of a cell cycle), we obtain after some manipulation 


(149) 

(150) 


E(„., r) - 

1 — eP'^ 


and the only term retained in the n —)■ oo limit is 


(151) 


TO/ \ n—>-oo 'f] 

]\L{m,r) - rriopt -^ 




(152) 


Different dynamic phases 

We are concerned with the extension of the generating function for the birth-death process over n cell divisions with the same 
rate parameters X,iy to the case where we have different dynamic phases described by parameters {Xi,i'i},{X 2 ,iy 2 }, ■■■■ The 
generating function for the birth-death process, without immigration, is simply hg from Eqn. |30[ We will write this expression 
in the form 


G{z,t\mo) 


/'Pz + QY° 

\Rz + s) ’ 


(153) 
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with coefficients 


P = 2^X{l + l'-2)-ri'{X + iy{l-2)) (154) 

Q = n'(X + iy{l-2))-2'^{Xl'+ iy{l-2)) (155) 

R = -Xri'(l-i) + 2^X{l + l' -2) (156) 

S = 2Xri'{l-l)-2^l{Xl'+ iy{l-2)). (157) 


If we now label these coefficients with an index r denoting the appropriate dynamic phase, so that, for example, Pr is Eqn. 
|154| with Xr, Uri Tir replacing A, n, we can write: 


hr = 


^overall — ^0 — 


Qr. 

Pj'h'p^\ ~\~ (^'p 
Rphpj^\ ~\~ Sp 

+ Qo 

Rqz H- Sq 


where Pp^ ...Sp are given by the recurrence equations 

Pr -- 
Qr - 

Rp - 
Sp -- 

with Pn = Pn^ ■■■! Sn = Sn- If we Write the matrix 


PrPr-\-l P QrPr-\-l 
QrQr-\-l P QrPr-\-l 

h^ph^P^'^ Sph^pj^"^ 

SpQp-\-\ P SpSp-\-\^ 


(158) 

(159) 

(160) 


(161) 

(162) 

(163) 

(164) 


M,. = 


Pr Qr 


Rr Sr 

the general solutions to these equations can compactly be given by 


(165) 


Po Qo 

Rq Sq 


11 “. 


(166) 




Birth-death-binomial extinction probabilities without balance and/or cell divisions 


The birth-death-binomial generating function is given by setting a = 0 in Eqn. 119 We set A = K-|-j^-|-ln2/T and only 
consider times immediately after cell divisions, hence setting t = 0 and I' = 1, giving 


_ / 2v{lk - 1) -I- T ^l^{z - 1 ){kt -\- VT(2lk - 1) -I- In2) -I- (1 -I- z - 2lkz){K + i/ + In2/T) 
’ \2L'{lk — 1) -1 T~^VJ:{z — l)(2Zfc — l)(Kr -I- -1- In2) -1- (1 -I- z — 2lkz){K -1-I- In2/T) 

where Ik = Extinction probability is thus 


(167) 


PsoOn = 0 ) = 


(Z(f — 1){kt -I- vT{2lk — 1) + In2) 


kt{ 211+^ vT{2lk - mi - 1) - In 2 - In 2 + In 4 


which reduces to 


(168) 


PBoim — 0) 


l-ll 


l + ll- 211+^ 


(169) 


for 1 / = 0, as given in the Main Text. 

In the absence of cell divisions, the birth-death generating function is simply Eqn. with a = 0. Setting z = 0 gives the 
general extinction probability 













23 




(170) 


Copy number balance can be enforced in the absence of cell divisions by taking the X = v limit, from which follows the 
generating function 


Gbd'm{z, t) 


z + vt — 

\ vt — vzt j 


from which straightforwardly follows the extinction probability 


PsD'Mi'm 



(171) 


(172) 


Other partitioning regimes 


We have derived results for the case where a binomially-distributed random number of agents is lost at each cell division. We 
now consider the case where this number is a fixed constant. We will denote this constant loss number by rj. In this case, 


Psirntyrriiy = 


(173) 

Mi 

[g{z, t)]"‘’'“ Ps{miymiy 


rrii^a—O 

yz,t)g{z,t)-'^g{z,tr^’^- 

(174) 

and so (j){z, t) = 

g{z,t)-^-, 

(175) 

S{g{,z,t)) = 

g{zA). 

(176) 


As 9, ^ and g take the same form as for the random loss case, the solutions for Zi and hi are the same as before. The difference 
(due to the different form of (jj) is in the second product in Eqn. which is now nr=i '^)~^- the Appendix we show that 
this factor takes the form of Eqn. 32 with Ai = vl{z — 1)(— 0 : 2 ) = B 2 = x"(i7 — Az), A 2 = \l{z — 1)(—X 2 )",Pyi = Xi,pb = 

(-X2),7 = -r]- 

The generating function in the case of deterministic subtractive inheritance involves (i) Eqn. |104 applied to the appropriate 


terms described in Eqn. 106 (ii) Eqn. 104 applied to Eqn. 113 as described; (iii) Eqn. ini and (iv) Hq from Eqn. |84[ and is 


Gn{z,t) 


- ix-*.)"+ - xz}} ( -"lyyiyr , 

I {xn{z -1) + H„-Xzn (s#Af;') „„ ) °^" 
(Ai'(z -1) +1(^ - A^)) j 

/ ly-X / GVyz - 1) + V - Xz X'^^ 

\AZ(z — 1) — Xz + u J \l^l'X{z — 1) + u — Xz J 


(177) 


Now we briefly explore two other inheritance regimes of potential biological applicability. In these cases we have not been able 
to obtain closed-form solutions for an arbitrarily large number of cell divisions: however, the appropriate recursion relations may 
be followed for as many divisions as required in order to obtain a closed-form solution for the generating function. 

First we consider deterministic partitioning of agents, where each daughter inherits exactly half of a parent’s population. In 
this case: 


Psim^yrmy = 


(178) 


[g(z,!)]"*'’“ Psirmyrmy 


rrii^a—O 

azMz,tr^'^/^-, 

(179) 

and so 4>{z, t) = 

1; 

(180) 

0{g{z,t)) = 

Vgiz,t), 

(181) 
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leading to the recurrence relations 


2;^ = y/g{zi-i,T) ; zi = ^g{z,t). (182) 

hi = g (^y/hi+i,T^ ;hn=g{z,t). (183) 

Next we consider the binomial inheritance of clusters of agents. We will assume that these clusters are of fixed size Uc. In 
this case, we consider the new variables Cf, = Ca = fni^al'^c (denoting the number of clusters before and after a cell 

division), and write 


OO '^i,b 


E E g'^'’‘^Ps{mi^a\'rni^b)Pi-i{mi^b,T\rrio) 

rrii^h—O rrii^a—O 
OO Cb / ^ \ 


Cb=0 Ca=0 

OO 


Cb=0 


= E(b + V) Pi-iincCb,T\mo) 


a 

Cb 


2 2 


= E + P^-l{mi,b,T\mo) 


rrii^b—O 


(184) 

(185) 

(186) 


The resultant generating function analysis yields a very similar outcome to that in Eqns. |15]20[ with the altered recurrence 
relation 


/I ^ 

^ I 2 2 



(187) 

(188) 


We have been unable to reduce the recurrence relations Eqns. 182p83| or Eqns. |187P88 to a closed-form solution for birth- 
death dynamics, but the corresponding problems may be solved for an arbitrary number of cell divisions by writing out the 
recurrence explicitly, thereby obtaining the generating function for a given number of cell divisions. The figure in the main text 
uses this approach. 
















